IMCI <- function(upper, lower, var.upper, var.lower, 
        conf.level=.95, eval.grid=.001){
    A <- (upper-lower)/sqrt(max(var.upper,var.lower))
    C <- seq(0,10,by=eval.grid)
    const <- abs(pnorm(C + A) - pnorm(-C) - conf.level)
    Cn <- C[const==min(const)]
    ci <- c(0,0)
    names(ci) <- c("lower","upper")
    ci <- c(lower - Cn*sqrt(var.lower), upper + Cn*sqrt(var.upper))
    if (is.na(ci[1])) ci <- c(NA,NA)
    list(ci=ci, conf.level=conf.level)
}
